function ro = ranktrans(rd, dim, opts)
% ranktrans  - return rank transform of data
%
% FORMAT:       ranks = ranktrans(data [, dim [, opts]])
%
% Input fields:
%
%       data        numeric data
%       dim         dimension along to transform (default: last)
%       ppts        optional settings
%       .meancenter flag, center around mean
%       .nozero     flag, do not "rank" zero samples
%
% Output fields:
%
%       ranks       1:size(data, dim) rank-transformed data
%
% Note: ties will get the mean rank.

% Version:  v0.7g
% Build:    9040106
% Date:     Apr-01 2009, 6:24 AM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~isnumeric(rd)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing first argument.' ...
    );
end
if ~isa(rd, 'double')
    rd = double(rd);
end
if nargin < 2 || ...
   ~isa(dim, 'double') || ...
    numel(dim) ~= 1 || ...
    isinf(dim) || ...
    isnan(dim) || ...
    dim < 1 || ...
    dim > ndims(rd)
    dim = ndims(rd);
end
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'meancenter') || ...
   ~islogical(opts.meancenter) || ...
    numel(opts.meancenter) ~= 1
    opts.meancenter = false;
end
if ~isfield(opts, 'nozero') || ...
   ~islogical(opts.nozero) || ...
    numel(opts.nozero) ~= 1
    opts.nozero = false;
end

% get elements that are Inf, NaN
bv = isinf(rd) | isnan(rd);
if opts.nozero
    bv = bv | (rd == 0);
end
mmm = minmaxmean(rd, 4);
rd(bv) = mmm(2) + 1;

% first pass, sort forwards
[rs, ri] = sort(rd, dim);
[ri, ro] = sort(ri, dim);

% to resolve ties (to mean rank) sort reversed array
sr = {':'};
sr = sr(ones(1, ndims(rd)));
sr = struct('type', '()', 'subs', {sr});
sr.subs{dim} = size(rd, dim):-1:1;
[rs, ri] = sort(subsref(rd, sr), dim);
[ri, rs] = sort(ri, dim);

% then build average
ro = 0.5 .* ro + 0.5 .* subsref(rs, sr);

% mean center?
if opts.meancenter
    rma = ones(1, ndims(ro));
    rma(dim) = size(ro, dim);
    ro = ro - repmat(sum(ro .* (~bv), dim) ./ sum(~bv, dim, 'double'), rma);
end

% remove bad values
if opts.nozero
    ro(ro == 0) = eps;
end
bv = bv | isinf(ro) | isnan(ro);
ro(bv) = 0;
